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Abstract We express the Crow-Kimura and Eigen models of quasispecies 
theory in a functional integral representation. We formulate the spin coher- 
ent state functional integrals using the Schwinger Boson method. In this 
formulation, we are able to deduce the long-time behavior of these models 
for arbitrary replication and degradation functions. We discuss the phase 
transitions that occur in these models as a function of mutation rate. We 
derive for these models the leading order corrections to the infinite genome 
length limit. 



1 Introduction 

The quasispecies models of Eigen |Hj and Crow-Kimura |3] are among the 
simplest that capture basic aspects of mutation and evolutionary selection in 
large, homogeneous populations of viruses. These models are a favorite entry 
point for physicists to evolutionary biology due to the phase transitions that 
they exhibit and their mathematical simplicity [Hj. 

While the models were originally defined in the continuous time limit, 
the first connections to statistical mechanics were made to the discrete-time 
versions of these models. In particular, the discrete-time Eigen model, which 
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entails the additional assumption of allowing only a single mutation at one 
point in time in addition to discretization of time, was shown to be equivalent 
to a particular type of Ising model |5j . A distinction between the bulk mag- 
netization and the observable, surface, magnetization was discovered in the 
discrete-time Eigen model l"> . and an analogy to the surface wetting phe- 
nomenon in condensed matter was made. Magnetization is a now-standard 
term in the physical quasispecies field for average composition. That is, the 
magnetization at site j on the genome is the average composition of the base 
at that site, averaged over all sequences in the population. A functional inte- 
gral representation of the discrete-time Eigen model was introduced through 
use of functional delta functions ^HI- With this representation, solution of 
this particular model was possible for viral replication rates that depend in 
an arbitrary way on distance from a single point in genome space. The closely 
related Crow-Kimura, or parallel or para-muse, continuous-time model was 
formulated as a quantum spin Hamiltonian in It was formulated as 

a functional integral and solved for viral replication rates that depend in 
an arbitrary way on distance from a single point in genome space in [111 
112113) . The discrete-time Eigen model in a sense interpolates between the 
continuous-time Eigen model and the Crow-Kimura model, because its limit 
as the time step becomes small is the Crow-Kimura model rather than the 
Eigen model. What distinguishes the continuous-time Eigen model from the 
other models is the possibility of multiple mutation events in an infinitesimal 
time step 

In this manuscript we seek to provide a detailed derivation for a functional 
integral representation of the continuous-time parallel and Eigen models of 
quasispecies theory. These representations are used to find exact solutions 
of these continuous-time quasispecies theories in the limit of large genomes. 
These results are used to exhibit the phase transitions that occur in these 
models as a function of mutation rate. The coherent states formalism that we 
introduce allows for the first time the expression of the full time-dependent 
probability distribution in sequence space of these quasispecies theories as a 
function of arbitrary initial and final conditions. In addition, we use the func- 
tional integral expressions to obtain for the first time the 0(1/ N) corrections 
to the mean replication rates in these models. Finally, we also use the func- 
tional integral expressions to find for the first time the 0(1/ VN) width of 
the virus populations around the most probable genomes in sequence space. 

The rest of the manuscript is organized as follows. In Section 2 we map the 
continuous-time parallel model onto a spin coherent state path integral using 
the Schwinger Boson method. We evaluate the theory for a general replication 
rate. We find the corrections to the infinite genome limit. In Section 3 we 
map the continuous-time Eigen model onto a functional integral, again using 
the Schwinger Boson method. While the functional integral appears more 
singular than in the parallel case, due to the presence of multiple mutations 
at a single time step in this model, we also solve this model for a general 
replication rate. We also find the corrections to the infinite genome limit. 
We discuss the width of the virus population in genome space in the large 
N limit and correlations in the field theory in Section 4. We also find the 
expression for the full time-dependent probability distribution as it depends 
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on initial and final conditions. We conclude in Section 5. Much of the detailed 
derivations are included in Appendices. 



2 Spin coherent state representation of the parallel model 

2.1 The parallel model 

In the parallel model, the probability distribution of viruses in the space of 
all possible viral genomes is considered. For simplicity, it is assumed that the 
genome can be written as a sequence of N binary digits, or spins: s l n = ±1, 
1 < n < N. Distances in the genome space are calculated by the Hamming 
measure: dij = (N — J2 n s n s n)/^- The probability for a virus to be in a given 
genome state, pi, 1 < i < 2 N , satisfies the parallel model differential equation 

= rtn - 2^ j=1 r oPi) + 2^ , : ''<»/', ■ w 

Here is the number of offspring per unit period of time, or replication rate, 
and /xy = [iA(dij — I) — N [iA(dij) is the mutation rate to move from sequence 
Si to sequence sj per unit period of time. Here A(n) is the Kronecker delta. 
The non-linear term in Eq. serves simply to enforce the conservation of 
probability, = 1- We can express the differential equation in a simpler, 

linear form 

with the transformation pi (t) = (t) / Y]^ qj (t) . The explicit form of the 
replication rate is = Nf(u), where u = (1/N) ^2 n s % n . 



2.2 The parallel model in operator form 



Motivated by the observation £Q that the parallel dynamics in Eq. J2J is 
equivalent to quantum dynamics in imaginary time, we express the model 
in an operator form. We define two kinds of creation and annihilation op- 
erators: a a (j), a£(j)) a = 1,2 and j = 1, ... ,7V. These operators obey the 
commutation relations 

[a a (i),ap(j)] = 







(3) 



These operators create either a spin-up state for a = 1 or a spin-down state 
for a = 2 at position j in the genome. While it might see more natural to 
introduce a single set of creation and annihilation operators to define whether 
the spin at position j is up or down, this approach leads to a non-Gaussian 
field theory even for a vanishing replication rate function. Use of two sets 



4 



of creation and annihilation operators leads to a Gaussian field theory, with 
non-quadratic terms stemming from the replication rate function. We find 
this second form of the theory more convenient for calculation. This second 
form, moreover, can be extended to the case where the sequence alphabet is 
larger than binary. Since the state is one and only one of the possible letters 
at position j, spin- up or spin-down in the binary alphabet case considered 
here, we will enforce the constraint that 



fi t i(j)]°[fito)] 1 |o,o> 



^2oi(j)a a (j) = 1 

a 

for all j. Thus, the state at site j is either 

|l ) 0) = [a t 1 (i)] 1 [4(i)]°|0,0) or|0,l> = [ 

Defining to be the power on a\(j) for spin state i, we can rewrite the 
parallel model dynamics as 

d 



(4) 



(5) 



dt 



P({n)})=n({n)})P({n}}) + » 



N r 

E 



(l-n))P({. 



■»[!•<{....„] 

We introduce the state vector 



■})-P({<}) 



•}) 



(6) 



$>(K})IK» 



8=1 

which satisfies the differential equation 



dt 



1=1 



dt 



(7) 



(8) 



We now write this Eq. (jSJ in operator form. First, we introduce vector no- 
tation for the creation and annihilation operators: a(j) = (oi(j),a2(j)) an( i 
8^(j) = (a|(j), a 2 (j)). Then we introduce operators 

Ti(j) =at( 3 > i ay) 
T = a^i) • a(i) 

with spin matrices 



1 

1 



CT2 



-1 

1 



CT3 



1 

-1 



(9) 
(10) 



Then the dynamics of Eq. JSJ) can be written as 

d 



dt 



m = -h\ii>) 



with 



H = Nf 



N 

E 



T 3 (j)/N 



N 

mE^c?) 



(ii) 



(12) 
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2.3 The field theoretic representation of the parallel model 



We convert this operator form of the parallel model into a functional integral 
by using coherent states. We define a spin coherent state by 



<J)) 



— f ,a t (i)- z O')- z *0')-a(j) 



l(o,o)Ci)> 



V n\m\ 

n,m— 

These coherent states satisfy a completeness relation 



3 = 1 



7T^ 



The overlap of coherent states satisfies 

(z'(j)|z(j)) = e ~2{ z "0')-[ z '0')- z 0')]-[ z '*0')- z *0')]- z 0')} 



(14) 



(15) 



Equation (jll|) enforces constraint |@J if the initial conditions obey the con- 
straint. To project arbitrary initial conditions onto this constraint, we use 
the operator 

N N 

3 = 1 3 = 1 

,2tt N 



n 

3 = 1 



^iJXi[^U)Mj)-i] 
2ir 



The probability to be in a given final state at time t is 

{no} 

Using the coherent states identity in a Trotter factorization, we find 



(16) 



(17) 



P({n},t)=E Jim <{n}| 



{no} 



n 

3 = 1 



dz* M (j)dz M {j) 



|{z M })({z M }|e 



-eH 



N 

n 



|{z M -l})({z M -l}| e ^ 



N 

n 



dzl(j)dzi(j) 



|{zi})({zi}|e 
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N 

n 

3=1 



dzo(j)dz Q (j) 



e -^P|{z })({zo}|{no})P({no}) 



JtaJ [Pz*Pz]({n}|{z M }> j £({z }|{n }}P({n„}) 

\{no} 



M 



Y[{{\z k }\e^ H \{z k ^})({ Zl }\e-^P\{zo}) (18) 



k=2 



For initial conditions that satisfy constraint |£fj 

({n}|{z M }) = Ue-' z " W ' ZM(j) n(j) • z M (j) 
3 

({z }|{n }) = n^ K0Tz ° Wz S(j) ' no(j) 



(19) 



Conversely, if the projection operator is needed for arbitrary initial condi- 
tions, we note 



P\Mf)) 



2tt N 

n 

27T W 



2tt 



r-Q^ e -a 3 - |(e a 3 - z)(i)) 

"'O .7=1 



(20) 



For initial conditions that satisfy constraint 10}, we may remove the projec- 
tion operator to find 

iV 

P({n},*) = lim / [Vz*Vz] £ P({n }) \{ n(j) ■ z M (j)z* (j) ■ n (j) 

{no} 3=1 



xe 2 

A/ 



{no} 



■lKO')-[«*C7)-«*-i(j)]-[(«);o")-«;-ia)]-«*-iO')} 



x ]^[ e 2i 

x ^ 2j=i WE (i) '! z t-i 0')- !] 
lim / [X>z*£>z] V P({n }) 



{no} 



x ]Jn(j) -z M (j) z*(j) • n (i) 

3=1 



-S[z*,zl 



(21) 
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where 



M 



5[z*,z] = z*(j) • z (j) + J^zJO') ' Mi) - Zfe-iO')] 



fc=i 

N 



£>£0>3Z fc -lOW 



M f 

fe=l k 

+^£[4(j>iZfc-i(.?)-l] 



+ 4/(zJ,z fc _i) 



(22) 



where e = t/M. To evaluate the expression involving A^/[{a^}, {a}], we use 
normal ordering. We define 



iV/[{at}, {a}] = N : /[{a*}, {a}] : +4/[{&t}, {a}] 



(23) 



where the notation : (•) : means in the operator expression for (•), place all 
of the {a^} to the left of the {a}. The additional terms that this operator 
commutation generates are collected in Af. We note that N : / : is O(N), 
whereas Af is 0(1). For example, for the quadratic replication rate f(u) = 



T^u 2 , we find 



tf/[^£ fit k>3 & tf)] = ^[^Eat(j)a 3 a(i)] 



3=1 



2 L iV 



3 = 1 



N 



2N 



a\i)a 3 (i)k(i)a\j)a 3 (j)a(j) 



i,3 = l 
N 



i,3 = l 



N 1 N N 

3=1 3=1 



N 



N:f(Z2&t(j)* 3 &(j)/N): + 

3 = 1 



(24) 



so that Af = 7/2 in this quadratic case. By induction, we can show that the 
general form of the commutation term is 



Af 



1 d 2 f(0 



2 

In the continuous limit, the probability at time t becomes 



(25) 



P({n},t) = J [Vz*Vz] £ P({n }) f[n 3 (t) ■ Zj (t) z*(0) • n,(0) e"^^ 

(26) 



{no} 



3 = 1 
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where we have switched the subscripts and arguments of the variables and 
where 



S[z*,z} = / dt' Vz'(t')[ft- Fl (j)+i(i)]z 3 (0+^ 

J ° 3 = 1 



N / dt'f 



N 



dt'Af[{2?{t')},W)}] 



(27) 



At long times, we find that e ^'Krio}) ~ e^ mt |n*) by the Frobenius- 
Perrone Theorem % , independent of initial conditions, where f m is the 
unique largest eigenvalue of H, and |n*) is the corresponding eigenvector. To 
evaluate this eigenvalue, we consider the matrix trace. We, furthermore, in- 
corporate the projection operator by the twisted boundary condition z (j) = 
e lXj ZM(j) arising from Eq. (J20J. We find 



Z = Tre~ tH P 

/•2tt 



N 

n 




lim / 

M->oo J 


3=1 







M N 

nn 

*!=ij=i 



dz* k (j)dz k (j) 



-S[z*,z] 



(28) 



where 



-S[z*,z] 



]J({z k }\e-^\{z k ^}) 



(29) 



fe=i 



with boundary condition Zo(j) — e lXj ZM(j)- The action is Eq. H22J) . without 
the initial Zq(j') • zo(j') term. Since the replication rate depends only on the 
total magnetization, the expression for Z can be simplified. In particular, we 
introduce = ^ YljLi z k(j) a 3 z k-i(j) to find, as discussed in Appendix A, 
that the partition function becomes 



Z = / [2>£Z>£]e- 



(30) 



where 



S[Z,$=N j det-f[W)]+£V)W) + ii-jjAf\-NhiQ (31) 



and 



= TrTe 



f* dt'[iJ.o-i+£(,t')o- 3 ] 



(32) 
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2.4 The large N limit of the parallel model is a saddle point 

The general expression of the parallel model partition function involves a 
functional integral. Using that N is large, this functional integral can be 
evaluated by the saddle point method. We impose the saddle point condition 
to find 



ss 
si 
ss 

~sl 



IrJ a^e-io 
= = N£ c - N— 



tide TrTeX)' ' W^i+L^] 

= = -iV/'(£ c ) + iVf c (33) 



£e,£c 

Evaluating the traces, we find 



Z<= tanht[ M 2 +g]^ (34) 



For large t we can solve Eq. (|34|) for £ c to find 



and evaluate Eq. I|9()|l to find 

lnQ~tJtf + g = t- 1 JL== (36) 
v 1 - Sc 

Using Eqs. i|35N3t)[l in Eq. I|31|) . we find that £ c is the value which maximizes 
InZ tig V Af 

= /KJ+^VH!-f + Y ( 37 ) 

This expression is the saddle point evaluation of the parallel model partition 
function. It is valid for arbitrary replication rate functions /. 

As an example, we calculate the error threshold for two different repli- 
cation rate functions. For our first example, we take the case of /(l) = A 
and / = otherwise. This case leads to the phase transition at A/fj, = 1. 
For A/fi > 1 a finite fraction p\ of the population is at £ c = 1, whereas for 
Aj \i < 1, all of the population is at £ c = 0. The fraction of the population 
at £ c = 1 is determined by the implicit equation pi/(l) + (1 — f>i)/(£ 7^ 
1) = h\Z/(tN)\^ , which gives p% = 1 — [a/A. For our second example, we 

consider the quadratic fitness /(£) = fc£ 2 /2 £Q. We find a phase transition 
at k/fi = 1, where the selected phase occurs for fc//i > 1 with an average 
magnetization given by £ c = ±\/l — (fi/k) 2 . The observable, surface magne- 
tization, u* , is given by the implicit expression /(u*) = hxZ/(tN)\^ so that 
= ±(1 — n/k). 
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2.5 0(1/ N) corrections to the parallel model 

We now evaluate the fluctuation corrections to this result. This procedure 
will determine the other 0(1/ N) contributions to the mean replication rate 
per site, (In Z)/(tN). We expand the action around the saddle point limit 



1 M 

s[€,Z] = s& s ,i D ) + - ]T 



d 2 s 



d 2 s 



k,l=l 



s & 5 & + 2 lOrw 



(38) 



We find 



d 2 S 



8 2 S 



d 2 S 



-Nef'(Q8 kl 
eNdu 



= -Ne 



Tr(T3e e ' i ~ fc ^ At0 ' 1+ ^ c<T3 ^(T3e e ^ M_ ' /_fc '^ A " :ri+ ^ c0 " 3 ^ 



+Ne 2 



Tr(73e £M ( Alcri+ £ c ''' 3 ) 

<-£ re eM(iJ,a- 1 +i a a 3 ) 



2 



(39) 



These terms, and the matrix trace, are evaluated in Appendix B. Solving the 
result of Eq. i|104|) for £ c in terms of £ c , we find 



InZ r—^ Af /"(&) 



1 fj, 



Using Eq. we find 
InZ 



i-f(tc)(i-e c ? /2 /» 



1/2 



tN 



= mc]+^vi^e c ^ f i 



1 M 



l-/"(€c)(l-e c 2 ) 3/2 /^ 



1/2- 



(40) 



(41) 



This is the expression of the parallel model partition function accurate to 
0(1/N 2 ). The expression is accurate for arbitrary smooth replication rate 
functions /. Shown in Figure □ is the comparison between this analytical 
result and a numerical calculation following the algorithm in pQ. 
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0.267, 



0.002 0.004 0.006 0.008 0.01 
1/N 



Fig. 1 The 0(1/N) shift in the free energy is shown (circles). Also shown is the 
prediction from Eq. 1)41^ (dashed line). We use f(m) = km 2 /2 with k/fj, = 2 and 
M = l. 



3 Spin coherent state representation of the Eigen model 

3.1 The Eigen model 

In the Eigen model, the probability distribution of viruses in the space of 
all possible viral genomes is considered, as in the parallel model. However, 
when a virus replication event occurs, the virus copies its genome, making 
mutations at a rate of 1 — q per base during the replication. The probability 
distribution in genome space satisfies 



dt 



2' 



3=1 



(42) 



Here the transition rates are given by Bij = q N ~ d ( l 'i) (\ — q) d ( l >i\ We define 
the parameter /i = A^(l — q)/q to characterize the per genome replication 
rate. We take /i = 0(1). We define q N = e~ Ml and note that in the large 
N limit, /ii — ■> [i. As with the parallel model, the non-linear terms simply 
enforce conservation of probability, and it suffices to consider the linear terms 
only 



dqj _ \ - 
dt ~ ^ 

3 = 1 



[B zj rj - '>,, Di_ <}j 



(43) 



As with the replication rate, the degradation rate is defined by Di — Nd(u) , 
where u = (1/^0X1 4- 
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3.2 The Eigen model in operator form 



Using the creation and annihilation operator formalism, the dynamics can 
be again written in the form of Eq. 1)11(1 , with 



H =U[qT a (j) + (l-q)T 1 (j)]Nf 



3 = 1 



N 



E T 3( fc )/^ 



Ds=l 



Nd 



N 



3=1 



N 



k=l 



Nd 



N 



3 = 1 



k=l 



Nd 



N 



3=1 



(44) 



where the last expression is valid for large N. Corrections to this expres- 
sion will come from \i\ ^ /x, the exponential not being exactly equal to the 
product, as well as normal ordering terms. We will address these corrections 
later. 



3.3 The field theoretic representation of the Eigen model 



We introduce the Schwinger spin coherent states. We consider the normal 
ordered form of the Hamiltonian, and first consider the expression : H :. We 
will consider the commutator terms later. We find P({n}, t) can be expressed 
as in Eq. I|21|l and the partition function Z can be expressed as in Eq. 128|) 
with 



s\ 



n<{z fe }k- £:H 1K-i}) 

fc=l 

M 

z*oti) ■ z o(J) + J2 Z IU) ■ i z k(j) - Zfc-iC?')] 



k=l 



k=l 
N 



N 



j2*tU)wk-i(j)/N 

3 = 1 



^z*(j)<7 3 zn(j)/iV 

3 = 1 



(45) 



For the partition function case, we have the boundary condition zo(j) = 

e ,;Aj z M (j)- We introduce £ k = ^ J2f=i z fc(j>3 z fc-i(j) and r) k = ^ J2f=i z I-(j>i z fc-i(j) 
to find, as discussed in Appendix C, that the partition function becomes 

Z = j \v£,V^][VfiVrf\e~ s ^ T '^ (46) 
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where 



S[t,£,fj,ri\=N I dt 
Jo 

-NlnQ 



'{- 



(47) 



3.4 The large N limit of the Eigen model is a saddle point 

The partition function of the Eigen model is represented as a functional 
integral. In the limit of large N, this integral can be evaluated by the saddle 
point method. The saddle point conditions are 



SS_ 
5£ 
6S 

SS_ 
6fj 

SS_ 

5r) 



= = Ni 



N 



TrTa^e 



J* dt'[r) c a 1 +£ c <j 3 ] 



TrTe-Zo dt '['' c,Tl+ ^ <J3 l 
Ne-^e^f{Q + Nd'(£ c ) + N£ c 



= Nr, c 



N 



TrTe 



f * dt' [fjc<Tl+tc 



= = -N^e-^e^"f^ c )+Nfj c 



Evaluating the traces, we find 



"21 1/2 



At long time, we find 



Vc 



— t a nht[f c +C 2 c } 1 / 2 



rf c +e c 



(48) 



(49) 



(50) 



from the second and fourth lines of Eq. I|49|l and from Eq. I|1U7|) . Using Eq. 
(I51JII in Eq. 1471) . we find that £ c is the value which maximizes 



In Z 

In 



= e-*e!W I= Gm)-d(tc) 



(51) 



This is the saddle point expression for the Eigen model partition function. It 
is valid for arbitrary replication rate functions / and degradation functions 
d. 
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As an example, we calculate the error threshold for two different replica- 
tion rate functions. For our first example, we take the case of /(l) = A and 
/ = 1 otherwise. This case leads to the phase transition at Ae~^ = 1 |14j . 
For Ae~^ > 1 a finite fraction p\ of the population is at £ c = 1, whereas for 
Ae - ^ < 1, all of the population is at £ c = 0. The fraction of the population 
at £ c = 1 is determined by the implicit equation pi/(l) + (1 — Pi)/(£ ^ 1) = 
\aZ/(tN)\^ , which gives pi = (Ae~^ — 1)/(A— 1). For our second example, 

we consider the quadratic fitness /(£) = l + fc£ 2 /2 [T2|. We find a phase tran- 
sition at fc//i = 1, where the selected phase occurs for k/fj, > 1 with magneti- 
zation £ c = ±V2[^l + ^ 2 (l + 2/£:)- 1 -[l 2 jkfl 2 j[L. The observable, surface 
magnetization, u*, is given by the implicit expression /(u*) = \nZ/{tN)\t , 
which for the Eigen model is a non-linear, transcendental equation even for 
the quadratic fitness case. 



3.5 0(1/N) corrections to the Eigen model 

There are 0(1/ N) corrections to Eq. (|51|l . The first comes from 



2N 



(52) 



The second comes from the normal ordering of / and d: 



(53) 



The third term comes from the approximation made in the last line of Eq. 

Ell. 



N N 

ii[ i +^c?)] = n :c ^ < 



: e 1 - 



(54) 



since : T a (j) n : = for n > 1 due to constraint The fourth term comes 



from normal orderingTi(j) inexp^^^ ^Ti(j)/N] andT 3 (j) in /E7=i T 3 (j)/N 
in /: 



N 



N 



j 



15 



di] 

N 



dm 



^ l Kn T,ij, : T:!j;T,ij, :]/N 



3 = 1 



+ 



de'^e^ 



N 



di] 



AT 



" /1 e^/fe) + |:e- M e^/'(a)K fc 



(55) 



where = i Ej=i z kO')°2 z fc— ■ Introducing this new field, we find 



k c = at the saddle point. Moreover, we find 



TrT<72e 



/ o t dt'[f; ( T 1 +C c 



= 



(56) 



TrTe 



The trace evaluates as (172) = 0. Thus this fourth, commutator term vanishes. 

There are also fluctuation corrections to Eq. (|51|1 . These are discussed in 
Appendix D. Using the results of Eqs. @5fr. ifKHjl. and (|129|l we find 



^ = e-^V^I m) - d&) + ±lf [l [1 - c'/&'] 1/2 



1 



(57) 



where the constants a', b' , and c' are defined in Eq. I|127|) . This is the ex- 
pression of the Eigen model partition function accurate to 0(1/N 2 ). The 
expression is accurate for arbitrary smooth replication rate functions / and 
degradation functions d. Shown in Figure [3 is the comparison between this 
analytical result and a numerical calculation following the algorithm in [Jj . 



4 The full probability distributions for continuous-time 
quasispecies theory 

In this section we derive the field theoretic expressions for the full probability 
distribution functions of the parallel and Eigen quasispecies theories. By 
use of the coherent states formalism, we are able to derive the distribution 
for arbitrary initial and final conditions. In the long time limit, the initial 
condition will not matter, as the system will reach a steady state. The final 
condition matters, though, because the weight assigned by the field theory to 
a given final condition is exactly equal to the probability that a given surface 
magnetization occurs in the population of viruses. 
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Fig. 2 The 0(1/N) shift in the free energy is shown (circles). Also shown is the 
prediction from Eq. 1571 (dashed line). We use f(m) — km 2 /2 + 1 and d(m) = 
with k/n = 2 and /i = 5. 



4.1 Field theoretic representation of the full probability distribution of the 
parallel model 

To obtain the full probability distribution P({n}, t), rather than simply the 
largest Frobenius-Perrone eigenvalue, for the parallel model from Eq. (|21[1 we 
add a term 



N 



(58) 



to the action. We find 



JV 



P({n},t)=£P({n })n 



{"ol 



dJ* 



(59) 



J*=J o =0 



Here a(j) = 2 — n(j) is the value of the spin at the final time, and ao(j) = 
2 — ^o(j) is the value of the spin at the initial time. Evaluation of this 
expression is carried out in Appendix E. Here we use the result of Appendix 
E for a couple different types of initial conditions. We define Qij to be the 

matrix element of the matrix Te-fo dt ^ r7l+ ^ t ^ a3 \ 

If, for example, we say that at t = 0, the spins are distributed randomly, 
but with a given initial surface magnetization, uq, then we take the terms 
in the multinomial expansion of (Qn + Q12 -*-•'<■ ' 1 - 1 A 
initial and final conditions: 



Q21 + Q22) that satisfy the 



P(u, t\uo) = lim 

M— >oo 



' M 

n 



ieNd£, k d£,k 
2^ 
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j 1 ^max[0, Af( ^ + " ) ] 

where ji + j 2 = iV(l + u )/2, ji + j 3 = N(l + u)/2, and ji +j 2 +j 3 +34 = N. 

Alternatively, if we take an initial condition with spin up and down equally 
likely ((u ) — 0) and define Q+ = Qn + Q21 and Q_ = Q 12 + Q22 we find 



P(u,t) = Um / 



n' ieNd^ k d^ k 
2ir 

.k=l 



4.2 The large N limit of the full probability distribution for the parallel 
model 

Since the full probability distribution is also expressed as a functional inte- 
gral, it can be evaluated by a saddle point in the large N limit. In the saddle 
point limit, this equation becomes 

fc=l 

1 ~T~ 11 1 — U 

+- ^-lnQ+ + — lnQ_ (62) 

with 

fit = /'(GO 

l+u itiZl [I + e^ci + e|/cr 3 ]cr 3 nfifc+i + ^l + e6 cr 3]ii+2i 

2 nf=i + e i i °i + 11+21 

1 - w Ilfri 1 [I + e^vi + e?icr 3 ]cr 3 ]li=fc+i [I + e V> a i + e6 cr 3]i2+22 
2 nfiil^ + ^l + e&cs] 12+22 

= -Z-<a 3 (fc)) + + — (trsCfc))- (63) 

In Appendix F we evaluate these expressions where the probability distribu- 
tion is large, in the Gaussian central region. 



4.3 The parallel model distribution function in the Gaussian central region 

Adding the terms from Appendix F together, we find that the probability 
distribution becomes Gaussian in the central region. That is Eq. I|62(l becomes 

InPM) 1 2 f f'(u,) \ 

— (const) — -ou — (64) 



N '2 V 2 /wt 
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Wc, therefore, conclude that 



2^u* 



<(»-<»» 2 > = w 



The fluctuation correction to the fitness is given by 

(/(«)) = /(u») + /'(«♦)<« - + \f"{uMH 2 ) +o(^ 



(65) 



(66) 



From expressions l|41f) and 1|65[) . we conclude that there is also a shift of the 
average magnetization for finite N: 



(u — u*) 



1 

WW) 



i - 



i-/"fe)(i-a 3/2 /M 



1/2' 



tiu*f"(u*y 
/'(«♦) . 



(67) 



4.4 An alternative derivation of the fluctuation corrections to the parallel 
model distribution function 



As an alternative, we may compute averages of the magnetization from the 
full functional integral expression for the partition function. Computing the 
surface magnetization From Eq. 161|l . we find 



(«)(*) 



£n' u'P(u',t) 
M 



= lim 



n 



Lfc=l 



2tt 



n 

Lfe=l 

x(Q+K~] + Q-[£]) 



lim 

A-f-KX 



2tt 



TV 



(68) 



This expression, however, is not easily calculated in the saddle point limit. 
Computing the fluctuation, we find 



(u 2 )(t) = 



M 



lim 

M— >oo 



n 



ieNd^d^ 



Lfe=i 

x(g+Kl + o- 



2tt 



(Q+j - Q-[fl) 2 + 4Q + [g]Q_[g]/iV 

(Q+[f] + Q-[e]) 2 
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/ lim 

/ M->oo 

This equation implies 



M 

n 

,fc=i 



ieNd^kdZk 



2tt 



(«-<«»' 




(69) 



(70) 



The term (1 — u 2 )/N is the variance of the spin at a given site 1 < i < N. 
The term ((5AQ/Q) 2 ), therefore, is N - 1 times the 0(1/N) correlations 
between spins at any two different sites 1 < i < j < N. Base compositions 
at different sites, while equivalent, are not uncorrelated. There is a 0(l/N) 
correlation between base compositions at different sites. 
The variance can be written as 



(u-(u)Y 



Jd(Su)e 



-Nf'(u»)Su 2 /(4u,)8S 



( § Q+-Q- \' 
\° Q+Q- J 



Jd(6u)e 



-N f'{u,)8u 2 /(iu») p 6S 



where 



SS/N 



HQ- 
1 + 



Q- 



"ml 



2 

— In Q + In Q- 



1 — u 1 — u 
— - — hi — - — 



(71) 



(72) 



We note that 



ss_ 

~N 



1 1 
2~ 



-f'{u*) 2 Su 2 + 0(5u 3 ) 



(73) 



Using Eqs. Q140JI . (|f4f|l . (fTT) . and lO we find exactly Eq. ijHiJ), which shows 
the consistency of the two calculations. This second calculation, however, 
makes explicit the contributions of spin-spin correlations to the fluctuation 
of the population in genome space. Shown in Figure 01 is the comparison 
between the analytical result of Eq. (|65|l and a numerical calculation following 
the algorithm in pQ. 



4.5 Correlation functions of the parallel model field theory 

Since the theory is Gaussian in the saddle point limit, we can evaluate cor- 
relations from the inverse of the Hessian in Eq. I|39fl . The inverse is given 

by 



-Nef"{£ c )I Nel 

eNI ~Ne 2 ^rrM + Ne 2 I 



(74) 
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Fig. 3 The 0(1/ N) shift in the variance of t he s urface magnetization is shown 
(circles). Also shown is the prediction from Eq. 1651 (dashed line). We use f(m) — 
km 2 /2 with k/fj, = 2 and n = 1. 



From Eq. (|1U(J|) we find 



1 /-/"(&) 1 



-l 



iV£ \ 1 b 4fc 2 e 2 +fc 2 + e / 

Inverting this matrix, and inverse Fourier transforming, we find 

I 1 1 ~ ^ c -2 £ |fc-i| A /b 2 ~/"(g c )M 2 /b 

^ \/i - (i - € 2 ) 3/2 /"(£c)//; 

{O&dO) = ^fci ^ Ofci 

N yRR!mw 

Ne N 

I /"(&) 1 ~ c -2e|fc-i| A /fc 2 -/"($ e ) M 2 /6 

* a/i - (i - C 2 ) 3/2 /"(6)/m ' 

(76) 



for k, I » 1 and fc, Z < M, where 6 = £t/\/l - £?. 



4.6 Field theoretic representation of the full probability distribution of the 
Eigen model 

The full probability distribution for the Eigen model is expressed as a func- 
tional integral in an analogous fashion to the parallel model. For the Eigen 
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model wc find 



P({n},t)= lim / 



' M 

n 

.fc=i 



ieNd£, k dt; k ieNdr) k drj k 



2tt 



2tt 



xe 



AT 



where 



x £P({no})n{[W)Uo>0')} 

{no} J=l 

M 

*;=i 

= f e / * dt ' [f ' ( *' )CTl+ ^*' ) ' T31 



(77) 



(78) 



Taking an initial condition with spin up and down equally likely, we find 



P(u, t) = lim 



' M 

n 

./c=i 



ieNd£ k d£ k ieNdr) k dr) k 



2tt 



2tt 



xe 
x 



TV 

JV(l + «)/2 



(79) 



4.7 The large iV limit of the full probability distribution function of the 
Eigen model 



Since the distribution function is expressed as a functional integral, it can 
be evaluated by the saddle point method. In the saddle point limit, the 
probability distribution function becomes 



In P(u,t) 
N 



M 

e^[e-^^/(a) - d(&) - - fj kV k} 
k=i 

l + ii, 1 + u 1 — u , 1 — 11 
-— ln — — ln — 

+ ^lnQ + + ^lnQ_ 



(80) 



with 



& = ^te(*0> + + ^%3(*)>- 



2 

1 + u 



%=/xe-" +Wfc /(a) 



>iW>- 



2 

1-u 



>i(*0>- 



(81) 
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Fig. 4 The 0(1/ N) shift in the variance of t he s urface magnetization is shown 
(circles). Also shown is the prediction from Eq. (1841 (dashed line). We use f(m) — 
km 2 /2 + 1 and d(m) — with k/jj, — 2 and jj, — 5. 



In Appendix G we evaluate these expressions where the probability distribu- 
tion is large, in the Gaussian central region. 



4.8 The Eigen model distribution function in the Gaussian central region 

The probability distribution function for the Eigen model becomes Gaussian 
in the large N limit. Using the results from Appendix G, we find that Eq. 
(|SU1) becomes 



In^) = _ 1 (f'M-fM 



(82) 



We, therefore, conclude that 



(u-(u)Y 



2/zm*/(u*) 



N[f'(u.) - d'(u*)] 



(83) 



From Eqs. (|57|l . and l|83|) we find the shift of the average magnetization 
to be 



1 



(U — M* > = 



Nf'(u*) 



b' 



1 - [1 - c'/b 



/d/2 



fj,u*f{u*)f"(u*) 
/'(tt*) — d'(u«) 



(84) 



Shown in Figure 0] is the comparison between this analytical result and a 
numerical calculation following the algorithm in jxj. 
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5 Conclusion 

We have derived exact functional integral representations of the Crow-Kimura 
and Eigen models of quasispccies theory. The functional integral representa- 
tion of these quasispecies theories is quite convenient because it allows us to 
obtain the exact infinite genome solution of these models as well as to obtain 
the finite length genome corrections. These exact results allow us to discuss 
the phase transitions that occur in these models as a function of mutation 
rate. The coherent states derivation of this functional integral also allows us 
to compute the full time-dependent probability distribution function of the 
population in sequence space, including the dependence on initial and final 
conditions. 

The functional integral representation leads to an interacting field theory 
for spin and mutation fields. In the limit of long genomes, the field theory 
becomes Gaussian. We have evaluated the theory in the infinite genome limit 
to give the mean replication rate for arbitrary replication and degradation 
functions. These results were used to exhibit the phase transitions that occur 
in quasispecies theory as a function of mutation rate. For smooth replication 
and degradation rate functions, we have evaluated the 0(1/ N) corrections to 
the mean replication rate. We have also derived the finite, 0(1/ VN) width 
of the virus population in genome space in this limit. 

The functional integral representation can be applied to generalizations of 
the quasispecies theories considered here. The extension of the present results 
to arbitrary replication and degradation functions that depend on distances 
from K points in the space of all possible genomes is straightforward with the 
introduction of overlap parameters |14|. With the Schwinger spin coherent 
state formalism, the extension to genomes with alphabets larger than binary 
0] is also straightforward, since the field theory in z* and z remains Gaussian 
due to constraint Q). 
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Appendix A 

To evaluate the partition function of the parallel model, we introduce £fc = 
jj ^2f—i z *kU) a 3 z k-i{j) with the representation 




(85) 



We thus find 



Z = 




Lfe=i 




Zfc(i)<S(i)fc<zi(j) 




= lim 

M— >oo 





(86) 
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The matrix S(J) has the form 



S(j) = 



(I 

-A 2 (j) I 
-A 3 (j) I 








(87) 



V ...-A M (j) I j 

where Ak{j) — I + efiai + effcO- 3 . We find 



det S(j) = det 



det 



M 

I~e^l[A k (j) 

fc=i 

EM - 
fc=1 u-i+e fc -3] 



Tr In 



Here the operator f indicates time ordering. The partition function becomes 

M 



Z = lim 

M— >oo 



n 



,fe=l 



2tt 



N 
3=1 



-Trlr 



e- iX *e 



= lim 

M— >oo 
JV 



n 

,fe=i 



2t? 



EM 



x n so*) 

3=1 



(89) 



where 



A/ 



QO') = Tr[][/ + e/i^ +e6a 3 ] 
fc=i 

= TwWa dt V CT i+5(*') CT 3] 



(90) 



Appendix B 

In this Appendix we evaluate the fluctuation corrections to the parallel 
model. The second term in the last derivative of Eq. (|39l) is given by 



Trcrqe eM ^ CTl+ ^ c<T3 ^ 



— tanht[ M 2 +Cl] 1/2 
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as t — > oo 



The first term in the last derivative is given for k ^ I by 
TrcT3e e ''~ fc " A " Tl+< ^' T3 ) ( 73e £ ( M- ''~' c '^ A1CTl+ 'f cCr3 ) 



(91) 



/ 



cosh(ti -* 2 )Vm 2 + Cc 



+ 4 2 + 4 2 cosh t vV + £? 



^ I ^ _ r (\tl-t 2 \-t)yf^+l 



+ 



» 2 +e c m 2 +4 2 



as t — > oo 



(92) 



where ti = e | Z — fc| and £2 = e(M — / — k\). For k = I, the first term in the 
last derivative is zero from the first line of Eq. (|90|l . We find 



where 



_£jL- = _Ar e 2 -^Mu + Ne 2 6 k 



(93) 



We let r^fc 
becomes 



M fei = e «[|M-2|i-fc||-M] > /^+^ (94) 
S^kyJ\eN) and 77^ = 5^k^/\eN), and the partition function 



Z = e - s (f<=.^) i im 



M . 



n 

Lfc=l 



idrjkdrjk 



2tt 



xe 



2 



(95) 



We let 77^, = iriky/f"(£ c ) and ^ = 7y fe / y/f"(£ c ) to obtain 



Z = e -5 ^'^ lim 



n 

,fe=i 



2tt 



xe 



2 St 



~ , >r- M kiilWi- e }"(£c) r l' k 7 l'i 5 ki 



Integrating over 77', we find 

M 



Z = £T s («-«<=) lim 



= e -S«c,«c) lim 

M— >oo 



n 

M 

n 



2tt 



2 12k 



Ski-* 11 J 2+ p } M kl +ef"(tc)Skl 



(96) 



U-l 

e -s(fo,fe) (detF) -i/2 



(97) 
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whcrc 



F kl = 6 kl - e ^fM Mkl + ef"{Q8 H 
= «w - e^M kl + ef"(£c)5 u 



(98) 



where a = /U 2 /"(£ c ) and 6 2 = y 2 + f 2 . We define F' , where F M = F' kl 
e/"(£ c )4z- We note that 

detF = det[J + e/"(&)J - e^M] 

= (det[J + e/"(£ c )I])(det[J - e±M[J + e/"^/]- 1 ) 
~(det[7 + e /"fe)/])(det[/-epM]) 
= (detF')e^ Hi+ef"m 
~ (detF')e Tr ^"^) / 
= (det 

To evaluate det F 1 we use Fourier space. We find 

4e 2 a 1 



(99) 



F'{k) = 1 



6 46 2 e 2 + P 



and 



detF' = JJ Ai = JJJ"(A) 



(100) 



(101) 



where k = 2m/M = 2nen/t, n = -M/2,...,(M - l)/2. In the limit of 
infinite M, we find 



det F' 



(sinht^/b 2 - a/b)' 



(sinh6t) z 



sinht 



1/2 



(sinhi^/i 2 + 

„ e _ 2 4V^l_^ +el ^ l2r (« c )/V^l] 1/2 ] 1/2 as t- oo (102) 

Using Eq. (|37|) . (|97|l . i|102|) . and (|99|) we find the mean replication rate at 
large time becomes 



1 M 



1 - 



w"(e c )(i-a 3/2 /v 



1/2' 



(103) 
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with 



[m 2 + £F 2 



tanht^ 2 + £;] 1/2 



(104) 



Appendix C 

To evaluate the partition function of the Eigen model, we define magne- 
tization and mutation fields with = J2jLi z %U) <7 3 z k-i (j) an d Vk = 

jj Ylj=i z fe0) cr i z fe-i U) to write the partition function for the Eigen model 
as 
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Af— >oo 



' Af 
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ieNd£ k d£ k 
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Af 
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ieNd£ k d£ k 
2ir 



' Af 
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2ir 



N 

: [VX]e~^l[[detS(j)r 1 



(105) 



The matrix 5(j) has the form of Eq. 1)8 7|) with = J + e^fcCi + eCfcC3- 

The determinant and constraint evaluate as in Eqs. I|88H89|I to give 



Z = lim 

Af— too 



' Af 

n 

,fe=i 



ieNdj k dj k 
2tt 



Af 
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ieNdrj k drj k 

2?T 



where 



,fc=i 

x e £jV SfcL 1 I 6 ^ 1 eMr,fc /«*)-<i(€h)-£*a-*7ki7*] 

JV 
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JV 



(106) 
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~ TrTe e ^fc=i [f ' fc ' T1+ ^ fc<:r31 
= TrTeX* dt '^ ( *' ),Tl+ ^*' ) ' T31 



(107) 
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Appendix D 



We here determine the fluctuation corrections to the mean excess replication 
rate per site, (In Z)/(tN), in the Eigen model. We expand the action around 
the saddle point limit 
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We find 
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2 / TrCTie eM ^ c<Tl+ ^°' 3 A / Tra- 3 e eM ( f '' :Cri+ ^ cCT3 ^ 



(1-40 



^Pj-gcA/CrjctTi+^cCra) 



The other two derivatives in Eq. (|108|l are zero. We find 



r p re eM(r )c CTi+4 c cr 3 ) 



tanht[^+g] 1/a 
as £ — > oo 



and 



TrCT3e e l' _A: l(' ,cCri+ ^ c0 ' 3 V3e e ( M_ l' _fe l^ r ' c ' T1+ ' 5c<J3 ) 



'J' re eM(T; c (Ti+^ c o'3) 



77 c 2 COSh(ti-f 2 )>/Tj2+£ 



?7 2 + C c 2 ?7 C 2 + coshtVWTH 



where ii = e\l — k\ and £2 = e(M 



fc|). We find 



d 2 s 



-Ne 2 _ n Vc ~ ML+Nf?6. 



m + e 



where 



ML = e 4\M-2\l-k\\-M\^j§+g 



The other traces evaluate as 

Trui e eM ^° CT1+ ^ c<T3 ^ 



'/<; 



'Pj.gCA/C^eCTi+^cCTa) 



[v 2 +e c ] i/2 
[v 2 +e c ] i/2 



tanh^+g] 1 / 2 
as t — > oo 



(109) 



(110) 



(HI) 



(112) 



(113) 



(114) 
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and 



Trcrie e l'" fc K , ' cCri+5c ' T3 Vie e(M_ l' _fc l) (r ' c ' Tl+ ^' T3) 
Tre eM ( f ' c<Tl+ ^ c ' T3 ) 



H COSh(i 1 -i 2 )0f+£2 



v 2 + e c v 2 + e c 2 cosh t vw+Tc 



n, 



2 , e c 



3 (|ti-* a |-t)Vfll+ef as t 



oo 



We also find 

, jy CTie t l(')cTl+f c CT 3 ) CT3e t2(rl c CT 1 +5 c CT3) r p rcr3 gt2(f; C 0'l+?cO-3) (Tie tl(r; c <Tl+CcO-3) 



(115) 



l^gtC^cri+lcCTs) 



r £ Te t(ri c a 1 +^ c a 3 ) 



V. 



,£ c cosh(ti - t 2 )^Vc +Cc 



^ + £ c 2 



fj 2 c + e c cosh*v¥+c! 



We, thus, have 



and 



3 (l*i-t a |-t)>/^ as t 



oo 



df] k df]i 



d 2 S 



= -Ne' 



m + il 



ML + Ne 2 6 k i 



8 2 S 



Ne 2 ^ c 



V 2 +C, 



(116) 



(117) 



(118) 



With these results, letting primes denote the fluctuation variables, and set- 
ting fii — /j,, the partition function becomes 



z = e -s&£c,vc,vc) i im 

Af^oo 



M 

n 

k=l 



ieNd? k d? k 
2tt 



M 

n 



2ir 



/c=i 



xe 



xe 



xe 



(119) 



We remove the primes, letting 77^ = rj k /\ftN, rj k = -ir] k /VeN, = ^ k /\ftN, 
and = —i£ k /\/eN, to find 



Z = e - S «c,£c,»7c,»7c) l im 

M— >cx> 



M 

n 

Lfe=l 



d£, k d£_ k 
2ir 



M 

n 



drj k drj k 
2ir 



xe 



>»(lc-l) f"c 



fe=l 
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xe 



xe 



(120) 



We let Xfc = (£,k,Vk)- We denote the summand that does not involve bars 



as -±x^B T .Bx fe ; that is B T B = e 



We let y fe = Bx fc and y k = B lT x fe . We note x fe • x fe = y fc • y fe . We let 







/-y _ J_ ( Vc ~ £c 



The partition function, thus, becomes 



Z = e -S(fe,£e,i?c,rjc) lim 



M , 

n dy k 
O.tt 



M ,_ 

2tt 



J Lfe=l 



"5 Z)fcll [2 ' 1 fc +J/2 ^ 2 ^ fc ' yfc ] 



M 

e -S(f c ,£ c ,77 c ,r, c ) Hm 



n dy k 
9ir 



,fc=l 



xe 



e -S(fc,6,77c,»7c) lim 
M— >oo 



n 

.k=i 

e -S(f c ,| c ,7jc,?7 c ) pll 



M 



d Yk 
2tt 



1 - T 771/ - 

"2 Z^fc,( = i y '= F fci y * 



(121) 



where 



Kl =S ki I-e 



ML 



kl 



: BC T CB T + eS kl BB T 



F'" + eS kl BB 



(122) 



We note 
det F" = (det 



/ + eBB T I - e 



M. 



kl 



ni + e c 



■BC T CB T 



(det[I + eBB T }) det 



- M 'kl_ j3r~ n n 



BC'CB 1 (7 + eSS 1 ) 



T\-l 



(detF'")[det(5 k i+eBB T 5 kl ) 



= (detF '") e Trln(«5 fci +eBB T 5 fci ) 

~ (detF"'^ 1 ™^ 
= (detF"')e msBT 
= (detOe tTrijTiJ 

= (det J F'") e teM<,,c_1) [/"fe)+A' 2 /(«.)]-w"(«e) 



(123) 
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To evaluate det F'" we again use Fourier space. We note 



where 



We find 



F"'(k) = 1- BC T CB T ^y- — 7y— 



A*c-"ew./(f c ) 



dctF'" =Y[dctF"'(k) 

k 

= TT dct (i - BC T CB T %, — 

11 V b' 4b ,2 e 2 + k 2 



n 



4eV 1 
6' Ab' 2 e 2 + k 2 



as i — > oo 



(124) 



(125) 



(126) 



where 



a' = TrSC 1 CB 1 = TrCB T BC T 



- d"(Qe^-^) - 2^/'(£ c ) +£m 2 /(£c) 
= b' 2 c' 



(l-e c )[f"^c)~e^^V^l d "^ c )] 



From Eq. fifty , we find 



(127) 



detF'" = 



waht^V -a'/V 
(sinh b'tf 

-2t[b'-[6' 2 -a'/6'] 1/2 l 



(128) 
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Using Eqs. (32), (52), (2D, (ED, (53), (S3), (EH, (123) , and fl2S) we find 
the mean excess replication rate at long times becomes 



InZ 

In 



= e^^/fc) - d(£ c ) - Uc ~ fjcVc + \/fj 2 c + e c 



N 



:b' 



1 - a'/b" 



1/2 



(129) 



Appendix E 

In this Appendix, we evaluate the functional integral expression of the full 
probability distribution function for the parallel model. Here Z has the form 



Z = lim 

M— »oo 



M 

n 

Lfe=l 

M W 

nn 

fc=0 j=l 



ieNd£ k d£ k 
dz* k (j)dz k (j) 



xe 



' M 

n 

,fe=i 



ieNd£ k d£ k 
2^ 



AT 

xe E" /"(^omOVW f[[detS(j)}- 



3 = 1 



where the matrix S(J) has the form 



S(j) = 

V o 

Equation l|59|) evaluates as 



( I 

-A 1 {j) I 
-A 2 (j) I 








-^Af(j) I J 



P({n},t) = Jim 



M 

n 



leNd^k 



2tt 



(130) 



(131) 



A' 



x ^P({n })n{[QO')]aoO>0-)} 

{no} 3=1 



eeE^Jiv/C^-AT^fc-A^+A/] 

(132) 
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where 

M 



Q(j) = n{I + ^oi +(£k<ra] 
fc=i 

~ 'f e c Z]fcli^ CTi +'f fc ' :r3 ] 

_ j> e / dt'[/iCTi+f(t')CT 3 ] (133) 



Appendix F 

In this Appendix, we evaluate the expressions necessary to determine the 
parallel model distribution function in the Gaussian central region. We define 
= limjv— >oo (w) and u = u* + <5ti. We note that = £(t') satisfies 

= u„i' = 

u,*'=* (134) 

We note that to <3((5w 2 ) 

+ = £(u», + (135) 

in the range 1<<', where 

St = — + fJ '~ f'( u *}/ u * Sy? (136) 
2/iti* (2^u*) 2 

because the differential equation 

^ = ,i(l + U )(a 2 (fc)) + + M (l - u){a 2 (k))- (137) 

is invariant under this shift, and the boundary condition £(i) = 77* + Su is 
satisfied to 0(Su 2 ) with the chosen value of St, Eq. (|136[) . We now evaluate 
Eq. (|n2Il to 0(Su 2 ). Using the shift property Ijl35|) and conditions Ijl34jl . we 
find to 0{Su 2 ) 

W) - f(Of(f) - = ( const ) - + 5t /M±^t±M 

+<^c/%)] - 5t u *nu*) + (u*+5u)fi(u* + 5u*) (i38) 



and 



1 + U 1+71 1 — U 1—71 1 + Ujf 1 + 1 — it* 1 — 77* 

— -z — In — — In — - — = — In — — In ■ 



2 2 2 2 2 

_ fem I±^_^! _J_ (139) 
l-Ti* 2 1-u v ; 
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Using the shift property (|135J) and boundary conditions l|134JI . we find to 
0(Su 2 ) 



[<5(u* + 5u)]- 



and 



[Q{u* + Su)]- 



-\*st 



1 + w* (1 - u* 1 + u* , . 

ot I — ^ — V + — ^ — / (u*. 



St 2 1 + 
~2 2"" 



(-2 /i « t /"K)+^ 2 + /'(»,) 2 ) 



(140) 



-A*5t 



1 — u» . 1 - , 

I — ^ — M / (w*) 



<ft 2 1 - 
"~2 2"" 



(2^/"K) + A1 2 + /'K) 2 ) 



(141) 



where A* = V7^+£c • 



Appendix G 

In this Appendix, we evaluate the expressions necessary to determine the 
Eigen model distribution function in the Gaussian central region. We note 
that to 0{St 2 ) 



£(u*+5u,t r ) =£(u*,t' + 5t) 
77(11* + 5u,t') = T](u*,t' + St) 



in the range 1 <§; t' , where 



St 



d'(u*) — f'(u*) + fiu*f(u*) + /iu 2 d'(u*) 



2/xu*/(tt*) 
because the differential equation 

dm 



U*[2/-iU*/(M*)p 



(142) 



<Su 2 (143) 



dt' 



(144) 



is invariant under this shift, and the boundary condition £(i) = + Su is 
satisfied to 0((5w 2 ) with the chosen value of St, Eq. (|143fl . 
Using the shift property (|142|) we find to 0(Su 2 ) 



= (const) - St[e-»+^f(£ c ) - d(&)] + St[f(u,) - d(u,)} 



StSu 



/'K)-rf'K) + / ./K)^ 

du 



(const) + 0(Su 3 



(145) 
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where 



dr\ 
du 



d 




*,t du 


J 



1 + u Q_ 1 — ui 



2 Q+ 2 Q_ 



d'(u*) - f'{u*) 



We also find 



-e(f)«o-^(* / M* / )] 



ii 



2 

dtdu 



(146) 



[/'(«*) - <?(«.) + u./"(u*) - u*d"(«*) + m/>*)] (147) 



We also find 



[Q(m« +fa)] a = e _ A . 5( 



1 + era* 



■ft 



1 — ctUi 



where 



di] 
IF 

d[ 
df 



-a^^l[f(u*)-d'(u*)} 



St 2 ( 1 — ra t c??y 
"2" V 2 ^ 

1 + OH 



1 + Qti, rf£ 

; 2 dt> 



[iff fa) 2 + (/'(«.) - d'(u*)f 



A* = £c£c + VcVc 



= — 2/ra*/(it*) 



/>♦)-*"(«♦) + /*/(«♦) 



(148) 



(149) 



We also have Eq. 1(135) . 



